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The density of an atom in a state of well-defined total angular momentum has a specific finite 
spherical harmonic content, without and with interactions. Approximate single-particle schemes, such 
as the Hartree, Hartree-Fock, and Local Density Approximations, generally violate this feature. We 
analyze, by means of perturbation theory, the degree of this violation and show that it is small. 
The correct symmetry of the density can be assured by a constrained-search formulation without 
significantly altering the calculated energies. We compare our procedure to the (different) common 
practice of spherically averaging the self-consistent potential. Kohn-Sham density functional theory 
with the exact exchange-correlation potential has the correct finite spherical harmonic content in its 
density; but the corresponding exact single particle potential and wavefunctions contain an infinite 
number of spherical harmonics. 



C/3 



X3 
O 

> 

O 

On 
OS 
-I— > 



I 

O 

o 



X 
S3 



I. INTRODUCTION 

Single-particle descriptions of electronic states and 
densities in atoms date back to their earliest models. 
Most of them involve the motion of individual electrons 
in some effective potential due the nucleus and the other 
electrons; Bohr's early analysis of some atomic spectra 
involved this idea With the advent of wave me- 

chanics, the idea took on the form of solving single- 
particle Schroedinger equations with this effective poten- 
tial. Prominent examples are the Hartree and Hartree- 
Fock (HF) approximations and density functional 
theory (DFT) Of these, only DFT provides in prin- 
ciple an exact description of electron densities with their 
proper truncated spherical harmonic content. In prac- 
tice, one is forced to adopt an approximate form for the 
exchange-correlation potential, such as the local density 
approximation (LDA). In carrying out such calculations, 
one computes the electronic states and effective potential 
iteratively, yielding a self-consistent potential and den- 
sity. 

The spherical symmetry of the nuclear potential yields 
states of well-defined angular momentum. However, ex- 
cept for S states, the resulting electron densities are gen- 
erally not spherically symmetric. (S'-states and their 
spherical densities present no problem and will not be 
further considered.) As we shall see below, in states of 
well-defined angular momentum quantum numbers L and 
Lz, the exact density n{r) may be decomposed in the fi- 
nite series 



n(r) 



E 



(1) 



where the n2/(r) are radial functions and l^jK^) 
spherical harmonics. 



The self-consistent densities obtained via the approx- 
imation schemes described above do not have this form. 
This form can be and often is obtained by introduction 
of a further approximation which, however, violates self- 
consistency: the effective potential may be spherically 
averaged, yielding single-particle states with good angu- 
lar momentum quantum numbers and a resulting den- 
sity of the form in Eq. |^. Such spherical averaging is of 
practical utility as it greatly reduces the numerical effort 
involved in carrying out the approximation schemes 
Nevertheless, the Hartree, HF, and LDA may all be ex- 
pressed in terms of variational principles, implying that 
the use of spherical averaging leads to an overestimation 
of atomic energy levels. To our knowledge, the quan- 
titative effect of this has only been checked in a small 
number of cases The effect is thought to be small 

since the resulting energies for many atoms are in quite 
good agreement with experiment 

In this paper, we will examine this, to our knowledge 
largely unexplored, issue in some detail. Using a pertur- 
bative approach, we will demonstrate that the inappro- 
priate spherical harmonic components appearing in the 
self-consistent density (without spherical averaging) are 
generally quite small. (An exact DFT calculation would, 
of course, yield the exact density with the correct spher- 
ical harmonic content.) We then develop a constrained 
search principle to modify the variational principles in- 
volved in the Hartree, HF, and DFT approximations to 
guarantee that the resulting density has the correct form 
of Eq. 1^. We show in the context of the Hartree approx- 
imation that this approach generates energies that are 
only slightly higher than those from the unconstrained 
approximation. 

It is interesting to consider in more detail the implica- 
tions of Eq. |l| for exact DFT. Being exact Q, it is un- 
necessary to introduce constaints to guarantee this "sym- 
metry" of the density. What is the angular symmetry of 
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the exact effective single particle potential entering the 
Kohn-Sham equations which guarantees that the density 
will have this form? A natural, but incorrect, guess would 
be that the Hartree and exchange-correlation potentials 
together sum to a potential that is spherically symmetric. 
In fact this is generally not true: the effective single parti- 
cle potential contains spherical harmonic components of 
all even orders. (A concrete example of this is presented 
in Appendix ^.) Indeed, it has been shown [|lO|,0 that 
a unique, spherically symmetric single-particle potential 
Voir) may be chosen to match the spherical component 
no(r) of the density; it is possible to formulate an al- 
ternative to the Hohenberg-Kohn theorem based solely 
on no{r) jl^. However, the resulting potential yields 
incorrect higher order components of the density n2i (r) 
(I > 0). Thus, there is a kind of complementarity: if one 
insists that the density have the correct truncated spher- 
ical harmonic content for an interacting state of well- 
defined angular momentum, the effective single particle 
potential will not have spherical symmetry and the sin- 
gle determinant model wavefunction will not have good 
angular momentum quantum numbers. By contrast, if 
one insists that the single-particle potential be spheri- 
cally symmetric, the correct spherical harmonic content 
of the density of the interacting state cannot not be re- 
produced. 

The remainder of this paper is organized as follows. 
In Section II, we first give the proof of Eq. 0. By a 
perturbative approach to the Hartree approximation, we 
demonstrate that the deviation of the density from its 
appropriate symmetry is actually quite small in a spe- 
cific example (a Helium triplet state), and comment on 
related results for the HF approximation and LDA. In 
Section HI, we formulate a constrained-search approach 
to single-particle approximations for the density, which 
we apply to the Hartree approximation and the LDA. 
We summarize our results in Section IV. Finally, two 
Appendices are included. In Appendix |^ we discuss a 
two-electron harmonic atom with interactions ("harmo- 
nium" ) and show explicitly that the density of its lowest 
triplet state cannot be reproduced by a non-interacting 
system in a spherically symmetric effective single-particle 
potential. Appendix ^ contains some details of the nu- 
merical calculations. 

II. SPHERICAL HARMONIC CONTENT OF THE 
DENSITY 

A. Finite Spherical Harmonic Content of the Density 

We begin by proving Eq. ^ of the Introduction by 
a standard application of the Wigncr-Eckart theorem. 
Consider an atom in a state \L,M > with total orbital 
angular momentum L and total azimuthal angular mo- 
mentum = M along the z-direction. We ignore spin- 



orbit coupling, so that the orbital and spin states of the 
atom may be specified separately. We are interested in 
the expectation value of the operator n(r) = 5{r — fi), 
where ri denotes the position of the ith electron. A use- 
ful decomposition of the delta function in this context is 

1=0 m=-l 

where fl represents an angular direction in spheri- 
cal coordinates. The set {Y^ifti), m = —I,— I + 
1, ...1} constitutes an irreducible tensor operator with 
respect to the angular momentum operator L, and 
obeys the Wigner-Eckart theorem ||l3|. It follows that 
(L,M|r;™(l1j)|L,M) is proportional to the Clebsch- 
Gordon coefhcient {LlMm\LlMm), which vanishes un- 
less m = 0, < I < 2L, and I is even. Substituting 
the expansion for the delta function into the expectation 
value of the density and using the above observation di- 
rectly yields Eq. |^. 

B. Infinite Spherical Harmonic Content in the 
Hartree and Hartree-Fock Approximations 

The decomposition of the physical density of an atomic 
state in spherical harmonics consists of a finite series. 
However, the Hartree, HF, and approximate DFT solu- 
tions do not produce densities with this property. For 
example, suppose we could find a finite decomposition 
for the density in the Hartree approximation, 

/— (/ even) 

The effective single particle potential contains a term of 
the form 

where A(0<A<l)isa parameter by which we may 
switch on the electron-electron interaction, which will be 
useful below. This term has a spherical harmonic decom- 
position with maximum I — Imax ■ The effective potential 
in the single particle Schroedinger equation multiplies a 
wavefunction (jj; for a single particle state of azimuthal 
quantum number m = this may be expanded as 

Cax 

H^=Y.yi'{r)Y,9{n). (2) 

l'=0 

When multiplied by the potential, the resulting products 
of spherical harmonics may be expressed as linear com- 
binations of single spherical harmonics Yj", with a max- 
imum non- vanishing contribution from I — Imax + ^'max- 



2 



The other terms in the Schroedinger equation however 
contain spherical harmonics of order I no greater than 
I'max- Thus, the Schroedinger equation cannot be solved 
by wavefunctions expressable in a finite spherical har- 
monic expansion (except for the trivial case of L' = 0). 

The density produced from these wavefunctions in gen- 
eral has no finite spherical harmonic expansion. One way 
to demonstrate this uses perturbation theory. The so- 
lution to the Hartree equations may be expressed as a 
power series in A; terms of higher order involve increas- 
ingly larger orders of spherical harmonics. When reorga- 
nized as a spherical harmonic expansion, all orders will 
occur with each coefficient a power series in A. It is not 
possible for these coefficients to vanish for arbitrary val- 
ues of A. 

As a concrete example, we analyze a two electron atom 
in a triplet spin state with total angular momentum 
L = 1, whose density is not spherically symmetric. Using 
perturbation theory in the electron-electron interaction, 
we compute the density in the Hartree approximation. 
The Hamiltonian for our system is 

ifl^2 ^™ \ri-r2\ 

where Z is the nuclear charge. The fully interacting sys- 
tem is given by A = 1, and we will formally develop our 
perturbation theory in powers of A. Alternatively, one 
may set A = 1 and consider an expansion of energy and 
density in powers of 1 /Z, which is equivalent to an expan- 
sion in A. Physically, one should thus think of the small 
A limit as the state of a highly ionized atom of large Z. 
The specific case we will focus on is = 2; thus A = 1 
and Z = 2 describes the helium atom. 

In the absence of interactions, the state of interest to 
us involves one electron in a Is state and one in a 2p 
state which we take to be in the m = state. It is easy 
to see that the density of this state satisfies Eq. 0: 

„(0)(.-') = |40)(r)p + |<)(f)|2 

= \RMr)Y„"m' + \R,,ir)Y,%n)f 

+cl,R,,{r)Y,"{n). (4) 

In Eq. ^ (/)q'^'' and 0^°^ are respectively the Is and 2p 
states, Rni{r) are hydrogenic radial functions, with n the 
principal quantum number and I the angular momentum; 
the coefficients c*j, are defined as 

cj, = j dCLY,°{CL)Y^{(l)Y^{n), 

so that Y^{n)Y^{Vl) = Y.^c]k'Y^{^)■ The coefiicients 
djf. are closely related to Gaunt coefficients commonly 



used in atomic structure calculations |9[], and have prop- 
erties similar to Clebsch-Gordon coefficients; in particu- 
lar, c^j. = unless i+j + k is even and \ j~k\ < i < \j + k\. 
It is these two properties that guarantee the density of 
the noninteracting state has the truncated form in Eq. 
|l|. The superscript, (0), in Eq. ^ denotes non- interacting 
quantities (A = 0). 

The Hartree 
approximation amounts to self-consistently finding two 
single-particle eigenstates 0o and 0i, of energies Sq, Si, 
for non-interacting electrons moving in an effective po- 
tential 

TT Ze^ o f r, , n(r') , , 

where the density is n(f) — |(/)o(r)P + |0i(r)p. To first 
order in perturbation theory, we may write that potential 
as 

11 J,. y \r-r'\ 

^_f^ + A;7«(r)+0(A2). 

ri 

The Schroedinger equation arising in the Hartree approx- 
imation may be solved within perturbation theory by ex- 
panding the effective potential, eigenstates, and eigenen- 
ergies in powers of A. The first order correction to the 
eigenstates satisfies the inhomogeneous differential equa- 
tion 

[Ho_er)]^W(f) = [e«_C/W]<^f(f), (6) 

where i = 0, 1, and the first order correction to the ener- 
gies are el"' ~ J d^r(j>f^* {r)U^^'' {r)(f)f^\f^. Because the 
density n^^^ (r) contains only even spherical harmonics, 
so will the potential U'-'^\ It immediately follows that 
the the wavefunction corrections (/)|^'' will have the same 
parity as the states from which they descend. Using 
the multiplicative properties of the Y^* gives 

=yo{r)Y„^{n)+y2{r)Y2°m, 

0« =2/i(r)rO(d)+2/3Mr3°(f^)> (7) 

where yi{r) are purely radial functions. 

Our development of the perturbation theory already 
illustrates one of the central points of this paper: we 
can see that the effective potential (which we have com- 
puted to first order in A) is not spherically symmetric 
and the wavefunctions arising in the Hartree equation do 
not have well-defined angular momentum. By expanding 
both sides of Eq. |^ in spherical harmonics and matching 
the coefficients for each I, the equations for the radial 
functions may all be written in the form 
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[ho{l) 



where 



ho{l) = - 



1 I d 2 d 
2m dr dr 



Mr) 



l{l + l) 
2mr^ 



(8) 



r 



and e 



(0) 



— £q°'' for even I, e^"' for odd I. The functions 
/; are easily computed, and the equations may be solved 
numerically. This calculation will be presented in the 
next section. Once the radial functions yi are obtained, 
the first order correction to the density in the Hartree 
approximation is found by adding the squared wavefunc- 
tions and collecting terms of order A. The resulting den- 
sity may be written in the form 



Spherical averaging is also a common practice in apply- 
ing the Hartree- Fock approximation and it is there- 
fore of interest to assess the extent to which Eq. |l| will be 
violated without such averaging. We again proceed per- 
turbatively. In addition to the direct potential C/*^^^(r), 
there is now a non-local exchange potential, to first or- 
der in perturbation theory, and the corrections to the 
wavefunctions take the form 



:il),HF _ ^HF 



i''^"^=yrir)Ynn). 



(10) 



with 



n^^\r) - 2cOoi?io(r)2;oM + 24,R2i{r)y,{r) 
4')(r) = 2cl^Rwir)y2{r) + 2c\^R2i{r)yi{r) + 2c?3i?2i(r 
n«(r) = 24i?2i(r)2;3M 

Note that to this order in A, only one "offending" 
spherical harmonic, Y^{yi), appears. However, all even 
spherical harmonics would appear in higher orders in per- 
turbation theory. Figure 1 illustrates the radial func- 
tions n'p {r) for the present model problem, as well as 
the analogous zeroth order densities nf^ (r) appearing 
in the spherical harmonic decomposition of the density 
for the non-interacting problem (cf. Eq. ^). Note that 
the magnitude of the offending spherical harmonic com- 
ponent is quite small (ratio of maximum contribution 
to root mean square density 3.86 xlO^'^), and that the 
densities n|^''(r) decrease very rapidly with increasing I. 
The reason for this is that at zeroth order, the density 
(of the non-interacting system) varies rather slowly as 
a function of the angular variable. When interactions 
are introduced, such a slowly varying potential has only 
a small amplitude for scattering electrons into high an- 
gular momentum states; the resulting density thus only 
has a small component of large I spherical harmonics. It 
is clear that this property is true at all orders in per- 
turbation theory: the effective potential entering at any 
order will always have a much larger Yq component than 
any other, leading to only small admixtures of high an- 
gular momenta in the wavefunctions. Finally, although 
we have illustrated this property in the specific context 
of a helium triplet P state, it should be quite general 
for atoms. Indeed, our model problem is in some sense a 
"worst-case" example; for larger atoms, particularly ones 
with many closed shells, the predominant spherical com- 
ponents of the density will be even larger. This helps to 
explain the success of using spherically averaged effective 
potentials in the Hartree approximation ||^. 



There is no Y§ term in the wavefunctions because there 
is a precise cancellation between the direct and exchange 
terms. The resulting density, remarkably, has precisely 
the right form - Eq. |l| - to first order, unlike in the 
Hartree approximation. In fact, one may demonstrate 
that the solution to the HE approximation reproduces 
the correction to the density exactly to first order in the 
electron-electron interaction. 
X / N Unfortunately, this good property of the solutions to 
the HE equations is limited to first order in A. This is 
(fiiost easily seen in the context of our model calculation 
for the He P state. The presence of a Y2 component 
to the density ensures that the effective potential seen 
by either electron has a similar component. For the (pQ 
state, this leads to a contribution proportional to Y^i^)^ 
as in Eq. |lO|. At second order in A, this necessarily pro- 
duces a component in the density proportional to Y^{Vl). 
However, the fact that the HE density has the correct 
form to order A indicates that the magnitude of the vio- 
lation will be even smaller than that found in the Hartree 
approximation. 



C. Harmonic Content of the Density in Density 
Functional Theory and Local Density Approximation 

The violations of Eq. |l| found in the Hartree and 
HE theories do not occur for the exact DET, which, by 
construction, produces exact densities. In practice, one 
must always introduce approximations for the exchange- 
correlation energy and potential. To illustrate the point, 
we consider a perturbative application of the local density 
approximation (LDA) to our helium P state example. 

The formalism closely parallels our perturbative ap- 
proach to the Hartree approximation. In LDA, we need 
to solve self-consistently a Schroedinger's equation with 
an effective potential given by V^^j + ^xc^^i where V^^y 
is given by Eq. | and 1/^^^ = SE^^'^lnir)]/ 5n{r). For 
the purpose of this illustration, we neglect the correlation 

1 /3 

contribution and take Vj^'J^^ir) w Vx{f) = -[f'^(?')] 
p^ . To first order in A, it is sufficient to replace n 
in Vj^'J^^ with n^"^. Because V^J^^ is not an analytic 
function of the density, it is important to recognize that 
when V^J^^ is expanded in terms of spherical harmonics 
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yj°(f2), the resulting series will involve all even values of 
I. This contrasts with the Hartree contribution, which at 
this order contained only I — and I = 2 components. 
In some sense this suggests LDA will lead to stronger 
violations of Eq. |^ than we encountered in the Hartree 
approximation. However, the effective potential we con- 
struct at first order in A is a slowly varying function of 
ri, so that contributions from large values of I to the 
wavefunctions are still quite small. 

When expanded in spherical harmonics, the correc- 
tions to the wavefunctions have a form very similar to 
Eq. 0, except (/ig^^ will now contain all even spherical 
harmonics, and (j)^^'' will contain all odd ones. Writing 

0-^^ = J2iy'2i+iY2i+ii'^), the equations satisfied by the 
y2i+i's are identical in form to Eq. ||, with a modified 
form for the inhomogeneous functions /; . Once the radial 
functions have been obtained, the first order correction to 
the density in LDA is given by n^^'^ (r) = J2i '^2/^ ('')^2/ (^) 
with 

41^ W = '2cl\2iRw{r)y2i{r) + 2c?^2;-ifi2i(r)2;2i-i(r) 

+2cl%+^R2iir)y2i+iir). 

In practice, the expansion of the density falls off so 
rapidly (see Fig. 3) with I that only the lowest few func- 
tions Un need to be computed. 

III. RESTORING THE SYMMETRY OF THE 
DENSITY 

A. Constrained Search Formulation 

The calculations in the above sections to some ex- 
tent explain why spherical averaging is successful in the 
Hartree and HF approximations, and in the LDA, when 
applied to atoms. Nevertheless, the averaging is basically 
ad hoc, and lacks a clear justification. From a formal 
point of view, spherical averaging has a dissatisfying as- 
pect: the minimization principles that are used to derive 
the three approximations are abandoned when it is intro- 
duced. Formally, a more consistent approach - especially 
for DFT - is to modify, or more precisely, constrain the 
wavefunctions searched in the minimizations in such a 
way that Eq. 1 is guaranteed. This has the advantage 
that the energies of the atomic states found will be lower 
than those found by spherical averaging. In practice, 
however, the energy lowering turns out to be quite small. 
Nevertheless, it is useful to explore constrained search 
methods for preserving symmetry properties of the den- 
sity because such violations are known to occur in other 
symmetry properties - particularly those involving spin 
- and may be responsible for more serious errors that 
arise in molecular calculations. The present formalism is 



a first example of how to consistently impose symmetry 
on an approximate single-particle scheme. 

As stated above, the Hartree, HF, and LDA equations 
are derived from minimization principles. In the Hartree 
approach, the energy functional is 

Em =< mo\^ > +y / d'rd'r'^^p^. (11) 

Here, IvE" > is a normalized wavefunction, Hq is the 
non-interacting electron Hamiltonian, and n{r) is the 
expectation value of the density in the state |\E' >. 
To generate the Hartree equations, one minimizes E[^] 
among orthonormal product wavefunctions 0. In den- 
sity functional theory, one adds an appropriate exchange- 
correlation energy E^c to the expression |ll|, and then 
searches for the minimum of the resulting energy 
||l5|,0. (This presumes the density is non- interacting 
u-representable, which we will assume for the states of 
interest.) For the exact exchange-correlation energy the 
resulting density satisfies Eq. |l|. Of course, the exact 
exchange-correlation energy is unknown, and, in prac- 
tice, one is forced to adopt approximations p^ . 

To constrain these searches to the subspace of states 
having a density of the form of Eq. |l|, we introduce a 
set of r-dependent Lagrange multipliers A2;(r). The con- 
straints that must be enforced are 

j dnn{r)Y^i{n) ^ 0, 21>2L, (12) 

where is L is the angular momentum of the state of in- 
terest. The function we need to minimize is 

E + E,, + j d\VR{r)n{r), (13) 

where 

l/^,(f)^ ^ A2,(r)r2l(f1). (14) 

1>2L 

After minimization of Eq. |l3|, we arrive at a single- 
particle Schroedinger equation 

[-^V^+Vs{^]U^ = eMr^. (15) 
For density functional theory, 

Vs{r) ^-^+e^ f dV'^^ + y.e([n(rO], rO + Vnif), 
r J \r — r\ 

(16) 

where V^c is the exchange-correlation potential. For an 
A^-electron atom, filling the lowest N eigenstates of Eq. 

leads to the density used in Eq. ^ so these equations 
must be solved self-consistently. The Lagrange parame- 
ters A2/ (r) of Eq. |l3| must be chosen to satisfy Eq. 
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One natural, but incorrect, guess would be that Vh 
simply removes the high spherical harmonics present in 
the other terms entering Vs , rendering a spherically sym- 
metric single particle potential. This is not possible ex- 
cept for the trivial case of S states. For example, in 
our model calculation of the helium P state, the lowest 
spherical harmonic component present in Vn is ^ = 4, 
which cannot remove the / = 2 component coming from 
the Hartree term. In fact, the single particle potential in 
general contains all orders of spherical harmonics. With 
the exact form of Vxc, Vr = 0, but Vxc itelf contains an 
infinite number of spherical harmonics. This is demon- 
strated in a specific soluble model (Appendix 0): two 
fermions in a harmonic trap, interacting via a repulsive 
quadratic potential. The exact eigenfunctions of this sys- 
tem may be written down explicitly, and in the Appendix 
we show (a) that the density of a P state cannot be pro- 
duced by a non-interacting electron system in any spher- 
ically symmetric potential, and (b) that the unique sin- 
gle particle potential reproducing this density (cf. the 
Hohenberg-Kohn theorem [^) in a non- interacting sys- 
tem contains all even orders of spherical harmonics. 



B. Perturbative Implementation 

Our proposed solution to the problem of producing 
densities that have an appropriate form for orbital angu- 
lar momentum eigenstates thus reduces to finding a self- 
consistent solution to Eqs. 12[ 1^, and The following 
is a practical procedure: (1) Obtain the self-consistent 
Kohn-Sham solution with the given i?2;c[n] and the self- 
consistent total potential Vs (r) . We expect that this will 
violate weakly the constraints with the offending 

density components n2L+2{'>'), n2L+4{r), . . . being small. 
(2) The restoring potential, VR.2L+2{r), Vr,2L+4 (?■),•■ ■ 
is determined by solving the equation 



IV. CONCLUSION 

This paper deals with the angular dependence of the 
electron density n(r) of an atom in a state of finite an- 
gular momentum L, both in the exact physical state and 
in various single particle descriptions. 

The spherical harmonic content of the physical den- 
sity, as a direct consequence of the Wigner-Eckart the- 
orem, is limited to even values of / < 2L. However in 
the Hartree, Hartree-Fock, and the various approximate 
forms of Kohn-Sham theory, the spherical harmonic con- 
tent of the density involves all Z-values (although com- 
ponents with / > 2L are small.) The exact Kohn-Sham 
effective single particle potential by definition reproduces 
the /-limited physical density; on the other hand, the po- 
tential involves all /-values. (This is documented for the 
case of an exactly soluble model of an atom with inter- 
acting electrons, the "harmonium" atom.) 

We show how the requirement, / < 2L, can be re- 
stored by a constrained search procedure using Lagrange 
parameter functions. Various numerical illustrations are 
presented. 

Somewhat analogous symmetry violations are known 
to arise in connection with the electronic spin quantum 
numbers. They may be susceptible to similar analysis 
and symmetry restoration. 

This work was supported by the NSF through Grant 
Nos. DMR 960452, DMR 9870681, and DMR 9976457, 
and by the Research Corporation. We thank Drs. D. 
Claugherty and Y. Meir for discussions of the finite spher- 
ical harmonic content of the physical density. 



-"2L-(-4(^) 
-"2L+6('') 



/ ^^.2L+2(^') \ 



\ 



(17) 



where K^^ is the submatrix (/, /' > 2L) of the linear den- 
sity response function Ki^ir{r,r') corresponding to Vs{r). 
(3) The new wavefunction and density, satisfying the con- 
straints (^2|) are determined from Vs{r) + Vii{r). 

Equivalcntly, this process can be carried out in terms of 
wavefunctions (see Appendix Perturbative energies 
for our helium P state example using the constrained 
Hartree approximation and LDA are presented in Table 
I along with comparable results for unconstrained and 
spherically averaged approaches. 



6 



APPENDIX A: DENSITY FOR A HARMONIC 
ATOM 

In this Appendix, we present a calculation of a P 
state for two spinless fermions trapped in a quadratic 
potential, interacting via a repulsive quadratic poten- 
tial. We call such a harmonic atom "harmonium" . This 
is not intended as a realistic model of a physical atom. 
But it shares with physical atoms their symmetry prop- 
erties and allows analytic calculations of wavefunctions 
and densities. Our main goals in this calculation are to 
demonstrate that (?) the exact density in this interacting 
state cannot be reproduced by non-interacting fermions 
in any spherically symmetric potential, and (ii) that the 
single particle potential that does reproduce the density 
contains spherical harmonics of all even orders. 

Our model Hamiltonian is 



H = - 



2m 



1 2 r 2 



(Al) 



Defining center of mass and relative coordinates fcm = 
(fi + ^2)72, f ~ (n. — ^2)72, this may be rewritten as 
a sum of commuting Hamiltonians, one for the center of 
mass coordinate [Hcm) and one for relative coordinates 



Hcm 
Hr 



2 

cm 



1 2 2 



(A2) 



where ^ = 2m, ujj^ ^ loq — The total angu- 

lar momentum operator may be written in the form 
L = LcM + Lfi, the sum of angular momentum operators 
for the center of mass and relative coordinates. Using the 
composition rules for angular momenta fl^ , it is easy to 
see that the P state of lowest energy is formed by putting 
the center of mass degree of freedom in an s state and the 
relative degree of freedom in a p state. Using the explicit 
forms for harmonic oscillator states, the wavefunction is 



1 n3/4 r T 

' exp 



2 

cm 



''CM 



2/2 -I 

"^'CM 



(A3) 



where Hiix) = 2a; is a Hermite polynomial, P — 
{fiujii)^^, and I'cf^j — {^,ujo)~^ . The density of this state 
is 

n{r) = [A + Bz^]e- 



{[A + isr^] + ^^lfr'Y,%n)}e-^'/^' (A4) 



^3/4 

B = (47r)3/2%Lc2 



(A5) 



The length scales appearing in the above two equations 
are given by L? = P + I'^f^j and = I'^j^^ + l^"^ , and 
C = [kIcmI]~^''^W^IV^- Note that Eq. ^ has the form 
required by Eq. |l|. 



A4 



is not derivable from 



We now demonstrate that Eq. 
a system of non-interacting fermions in a spherically sym- 
metric external potential. To show this, suppose the den- 
sity was derivable from such a potential. Then the two 
occupied single particle states would necessarily have the 
form 00 (f) = s\r)Y§[n), 0i(f) = P{r)Y^[n). The sum 
of the squares of these gives the density; matching this 
to Eq. A4 gives explicit expressions for S{r) and P(r), 



S{r) = Va^c-'' 
Pir) 



- B 



(A6) 



It is interesting to notice that the functional forms of 
S and P are perfectly compatible with a state of non- 
interacting fermions in a harmonic trapping potential. 
However, the normalizations of S and P are not correct. 
For example, S{r) is properly normalized if and only if 
A = {'kL'^Y_[\ An examination of the explicit expression 



for A, Eq. A5, reveals that this is the case only if wi = 0; 
i.e., the repulsion vanishes. Thus, for non-interacting 
fermions, no spherical potential will reproduce the inter- 
acting density. 

On the other hand, a non-interacting potential that 
is not spherically symmetric can be found to make the 
density of two non-interacting fermions take the form of 
Eq. A4. We again demonstrate this by using perturba- 
tion theory. Because the density is axially symmetric - 
i.e., it may be written as a function of r and z ~ we look 
for an effective potential which is also axially symmetric. 
In the body of this work we have essentially expressed 
the z-dependence of densities and potentials in terms of 
spherical harmonics. However, in this Appendix because 
harmonic oscillator wavefunctions have a number of use- 
ful algebraic properties, we expand instead in powers of 
z. 



The form of Eq. A4 suggests that the effective single- 
particle potential that reproduces the density has the 
form 



■6V{z), 



with 



where We// — 1/ ^iL^. We will perform our perturbation 
theory around a Vpff with 5V = 0; this is slightly dif- 
ferent than working around the non-interacting state, as 
the length scale L is modified by interactions. Neverthe- 
less, it is clear that 5V must be small if the interaction 
strength is weak. For this form of the potential, it is also 
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clear that the two single-particle states must have the 
form 



n even 



(A7) 



with i = 0, 1, and ipo the ground state of a one- 
dimensional harmonic oscillator with frequency ujef / , and 
for SV — 0, Xoiz) = ipo{z) may be taken as a harmonic 
osciUator ground state and Xii^) = "01 (-z) as the first 
excited state. If the electrons were non-interacting, we 



would necessarily have in Eq. A4 



A A r 1 n3/2 



(A8) 



It is convenient to parameterize the effect of interactions 
in terms of the deviations of A and B from these values, 

6n{f) = [SA + SBz^]e-'^'/^' 

with 6 A = A — Aq, 6B = B — Bq. Denoting the first order 
corrections to Xi(z) as Sxi{z), it is easy to demonstrate 
to lowest order in perturbation theory that 

[5A + 5Sz2]e--'/2^' - Sxoiz) + V2ySxiiz), (A9) 

where {SA,SB) = {nL^f^'^{SA,6B)/2. With some alge- 
bra Eq. may be written as 



6xo{z) + V2^Sxi{z) ^ PM^) 



(AlO) 



where /3 = ^{ttLY^H^BL^ - M], and i>2{z) is the 
n — 2 harmonic oscillator state. 

We now expand 8xq(z)^ (5xi(-z) in harmonic oscillator 
states: 



^Xq{z) = ^ C2-ai)2-a{z) 
oc 

^X\{z) = ^ C2„+l-02«+l(2:)- 



ri=0 



The coefficients c„ obey the recursion relation 
\Jn^ lc„+i = P8n.2 - c„ - V"^c„_i 



(All) 



(A12) 



for n > 2. In addition, cq, ci = since (5xo(-2), (^Xil^^) 
must be orthogonal to XO) Xi- 

Perturbation theory defines (5xoi (^Xil^) ™ terms of 
A useful expansion for SV is 



bV{z)= J2 



n even 



(A13) 



First order perturbation theory yields a linear relation 
between the Cn's and the u„'s. 



{n - l)uJeff 



n odd 



(A14) 



Eqs. 



A12 



and 



A14 may be combined to form a set of 
Any choice of V2 will generate 



recursion relations for v„ 
an entire set of w„'s, whose magnitude in general grows 
rapidly with n. The resulting 5V{z) is ill-defined. How- 
ever, for a given (3 there is a unique choice of f 2 for which 
\vn \ uniformly decreases with increasing even n. (Note w„ 
vanishes for odd values of n and V2m alternates in sign.) 
This choice may be found as follows: select a large value 
of N and set w„ = for n > N . The recursion relations 
for w„ (0 < 71 < N) then may be written as a matrix 
equation which may be solved numerically with little dif- 
ficulty. Fig. P illustrates Vn for /? = 0.1 and several 
increasing values of A^. As may be seen, the w„ converge 
to a unique sequence as TV —> 00, and the limiting values 
vanish rapidly with increasing n. 

The potential we have found is the unique single par- 
ticle potential that produces the density of the harmonic 
atom model, to first order in perturbation theory, as re- 
quired by the Hohenberg-Kohn theorem ||^. It involves 
all even powers of z, or, equivalently, all spherical har- 
monics Y^i- This demonstrates that, to reproduce the 
physical density, involving I = ane I = 2 only, requires, 
in the absence of interactions, an external potential in- 
volving all even values. 



APPENDIX B: PERTURBATIVE CALCULATION 
OF DENSITIES AND ENERGIES FOR 
CONSTRAINED ENERGY FUNCTIONAL 

In this Appendix, we discuss a concrete example of the 
ideas developed in Section III, providing details of how, 
for the helium P state, our perturbative Hartree approx- 
imation is modified by the introduction of the constraint. 

We begin with the single particle Schroedinger equa- 
tion, Eq. nsl with single particle potential 



Vs{r) 



Xe' 



, n{r') 



+ \VR{r) (Bl) 



To lowest non-trivial order in A, the corrections to the 
wavefunctions (jJf'^ obey Eq. |^. The form of t/^^-* must 
be modified to include the restoring potential, 



t/(i)(f)_{/)^^^(f) 



Note that to this order in perturbation theory, the I = 
and I = 2 components in a spherical harmonic expansion 
of u'^^ (r) come only from the Hartree potenti al, a nd so 
are identical to the ones encountered in Section [IB. The 
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higher I components of u'^^ come only from Vr. U^j^\r) constrained Hartree approximation. As in that approxi- 
thus can be expanded as mation, the introduction of the constraint introduces ht- 



tie change in the energy. 



+V4{r)Y^in)+veir)Ye°in) 



The radial functions are determined fully by the 

nonintcracting electron density, and the functions u„(r') 
come from Vr. 

Because all even values of / appear in the decomposi- 
tion of , the corrections to the wavefunctions contain 
all spherical harmonics: 

0^^^ ^ yoir)Y,'in) + y2{r)Y^"m + yi{r)Y,°m + - 

0W = yi(r)y°(f1) + y3{r)Y,'m + y5{r)Y,"m + ... . (B2) 



Noting that the radial functions yi (r) fall off very rapidly 
with we retain only / < 4 in the calculations that follow. 
The radial functions yi obey Eq. ^, with 



Mr 
fi{r 
f2{r 
h{r 



ir)Rio{r) 
cioV4:{r)Rio{r). 



(B3) 



Only /a and are modified by Vr; it follows that the 
2/0, yi, and y2 are identical to the results of the Hartree 
approximation without the constraint. The equations to 
be solved are closed by requiring the highest spherical 
harmonic component retained, Z = 4, to vanish to first 
order in A in the density: 

n«(r) = 4iRio{r)yi{r) + ctMr)y3{r) = 0. (B4) 

The constrained Hartree state is found by solving Eqs. 
^, B3, and B4. The radial functions ynir) resulting are 
presented in Fig. 2. We present in Table 1 the energies 
obtained by the Hartree approximation, the spherically- 
averaged Hartree approximation, and the constrained 
Hartree approximation. As can be seen, the differ- 
ences among the three approaches only arise at order 
and are quite small. Nevertheless, the con- 
strained Hartree approximation yields an energy consid- 
erably closer to the unconstrained Hartree approxima- 
tion than the spherically averaged one, indicating that 
imposing the symmetry by spherical averaging raises the 
energy considerably more than necessary. 

Finally, for comparison we also present in Table 1 anal- 
ogous energies for the results of a constrained, perturba- 
tive LDA calculation. The method for computing the 
wavefunctions is, mutatis mutandis, the same as for the 
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TABLE I. Zeroth {E^°^), first (AS'^^), and second (A^S'^^) 
order contributions to the He 2P state energy calculated 
in perturbation theory in the electron-electron interaction 
strength A by various methods. i/=Hartree approximation, 
H—Sph= spherically-averaged Hartree approximation, H—C 
— constrained Hartree approximation, LDA = local density 
approximation, LDA — C= constrained local density approx- 
imation. Energy units are e^/aB with os = fi^/me^ the hy- 
drogenic Bohr radius. 





H 


H-Sph 


H-C 


LDA 


LDA-C 


SCO) 


-2.5000 
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-2.5000 


-2.5000 


-2.5000 




1.3063 


1.3063 


1.3063 


0.5439 


0.5439 




-0.4059 


-0.4039 
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-0.1140 


-0.1139 
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FIG. 1. Radial functions of the spherical harmonic decom- 
position of the density in the Hartree approximation in a 
helium P state. Zeroth and first order corrections (desig- 
nated by superscripts) in the electron-electron interaction arc 
illustrated. Absolute magnitudes fall quickly with increasing 
spherical harmonic index / (designated by subscripts), (a) 
4°)(r), n«(r); (b)nf (r), n«(r); (c) n^^\r). 
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FIG. 2. Radial functions yi{r) produced by the constrained 
Hartree approximation to first order in perturbation theory 
for a helium triplet P state (see text), (a) yo{r), yaij); (b) 
y2{r), ysir), 2/4 (r). Note that the Z = 0, 1, 2 contributions are 
identical to the results of the standard Hartree approximation; 
I = 3,4 results are modified by the constraint. 




FIG. 3. Approximate coefficients v„ in expansion of 5V 
given by Eq. A13 for /3 = 0.1 (see text). v„ is assumed to 
vanish for n > N with different values of A'^ given in the figure. 
The expansion coefficients converge to a unique set of values 
as N —f oo. 
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